A free boundary (or moving boundary) generally exists between fluid and structure, and helped by the improved performance of recent computers and improved numerical simulation techniques such as the finite element method, the fluid structure interaction simulation is used not only in fundamental science but also in applied fields such as biomechanics. Hence, it has become possible to understand phenomena under more realistic problem settings. The simulation with respect to the fluid structure interaction problem often uses the IT (Interface-Tracking) type ALE (Arbitrary Lagrangian-Eulerian) method, as described in Nobuyoshi Tosaka, Hideaki Miyata, Yoji Shimazaki, Takashi Nomura, Masayuki Shimura, and Katsumori Hatanaka, Numerical Fluid Mechanics Editing Meeting Edition, “Numerical Fluid Mechanics Series 4 Moving Boundary Flow Analysis”, University of Tokyo Publisher, Tokyo, 1995 (hereinafter Tosaka et al.) and Q. Zhang and T. Hisada, “Analysis of Fluid-Structure Interaction Problems with Structural Buckling and Large Domain Changes by ALE Finite Element Method”, Comput. Methods Appl. Mech. Engrg., Vol. 190, pp. 6341-6357, 2001 (hereinafter Zhang et al.), for example.
FIG. 1 is a diagram for explaining a numerical simulation using the IT type ALE method. An algorithm or program according to the ALE method is executed by a computer system 1 illustrated in FIG. 1. The computer system 1 includes an input device 2, a computer 3 (or a processor) 3, and an output device 4. The input device 2 inputs to the computer 3 information including physical constants of the fluid and structure (for example, viscosity coefficient, density, Young's modulus, etc.) and initial condition, computing region (or domain), and mesh information. The output device 4 outputs information computed in the computer 3, including numerical information on the degree of freedom of the fluid pressure at each discrete time step, changes in the velocity or structure, visualization information, and the like.
The algorithm or program according to the ALE method includes a mesh information read part 6 to read mesh information from the input device 2, a graph information forming part 7 to form graph information called connectivity among nodes of discretised computing regions (that is, a node list forming finite elements), an initial condition set part 8 to set the initial condition, and a main time development loop part 9 to repeat N time steps.
FIG. 2 is a diagram for explaining an example of the structure of the main time development loop part 9 according to the IT type ALE method. The main time development loop part 9 includes a FSI part 11 and a mesh control part 12. The FSI part 11 forms the unit that simulates physical phenomena. The FSI part 11 includes a matrix creating routine 111 to create a small matrix and an equivalent nodal force of each element of the structure in correspondence with the time discretization technique, a matrix creating routine 112 to create a small matrix and an equivalent nodal force of each element of the fluid in correspondence with the time discretization technique, and a boundary condition applying routine 113 to apply the boundary condition of the structure and the fluid. The FSI part 11 further includes a matrix assembler routine 114 to create a matrix equation with respect to all degrees of freedom by restructuring the contribution from the small matrices of each element obtained from the routines 111 through 113, a linear solver 115 to solve and evaluate a variation of dependent variables (pressure, velocity, displacement, etc.) of the fluid and the structure every time a Newton-Raphson loop within the FSI part 11 is executed, and an output part 116 to output, by displaying, for example, the variation of the dependent variables and the evaluation thereof obtained by the linear solver 115.
If the displacement of the structure is denoted by “d”, |d|≠0, and the node mismatch occurs in the boundary of the fluid and the structure every time the Newton-Raphson loop is executed within the FSI part 11. The mesh control part 12 within the main time development loop 9 has the role of recovering from mesh mismatch (or mesh inconsistency) caused by the node mismatch, and forms an auxiliary part to enable the numerical solving. More particularly, the mesh control part 12 recovers from the mesh mismatch state to the mesh matched state by using the displacement of the structure as the boundary condition and treating the mesh of the fluid region as an elastic body. For this reason, the mesh control part 12 includes a matrix creating routine 121 to create the small matrix for each element of the structure, a boundary condition applying routine 122 to apply the boundary condition of the displacement of the structure, a matrix assembler routine 124 to create the matrix equation with respect to all degrees of freedom, and a solver 125 to evaluate the state of the mesh matching (or consistency of the meshes, hereinafter referred to as “mesh consistency”) every time a loop is executed within the mesh control part 12. If the target of the mesh control simulation is treated as a non-linear elastic body, the loop of Newton-Raphson method is executed iteratively in the mesh control part 12. In contrast, if the target of the mesh control simulation can be treated as a linear elastic body, the process escapes from the loop when the Newton-Raphson loop is executed once within the mesh control part 12.
According to the ALE method, the mesh control part 11 must carry out the mesh control described above, however, the mesh control has a possibility of failing particularly in a large-displacement problem. When the mesh control fails, the execution of the simulation stops, and an approximate solution to be obtained may not be obtained. In addition, in a problem in which the structures contact each other, the meshes will definitely fail (because a fluid mesh having zero volume is generated by the contact), and the application of the mesh control is limited. Moreover, because the amount of computations to be executed by the mesh control part 11 is relatively large, a load of the computations on the mesh control part 11 (that is, the computer 3) is relatively large. Accordingly, the conventional fluid structure interaction simulation method using the IT type ALE method is dependent on the mesh consistency between the fluid and the structure, and may not treat the large-deformation problem and the problem in which the structures contact each other. Furthermore, the load of the computations on the computer is relatively large because of the relatively large amount of computations required by the mesh control.
On the other hand, methods of performing the fluid structure interaction simulation using the IT type ALE method without carrying out the mesh control have been proposed in Japanese Laid-Open Patent Publications No. 9-245080 and No. 2000-352545, Toshiaki Hisada and Takumi Washio, “Numerical Study Related to Fluid Structure Interaction Simulation Method of Aortic Valve”, The Japan Society for Industrial and Applied Mathematics, Vol. 16, No. 2, pp. 36-50, 2006 (hereinafter Hisada et al.), and J. de Hart, G. W. M. Peters, P. J. G. Schreurs, F. P. T. Baaijens, “A three-dimensional computational analysis of fluid-structure interaction in the aortic valve”, J. Biomech., Vol. 36, No. 1, 103-112, 2003 (hereinafter de Hart et al., for example. However, these proposed methods are of the IT type that depends on the mesh consistency, and the meshes of the fluid become relatively small if the structure has a complex geometrical shape, for example. Hence, the amount of computations associated with the small meshes becomes relatively large in the case of the structure having the complex geometrical shape, and the load of the computations on the computer is relatively large.
When treating the large-deformation problem or the problem in which the structures contact each other, or when reducing the amount of computations, it may be desirable to perform a fluid structure interaction simulation that does not depend on the mesh consistency. However, the pressure and the velocity gradient change discontinuously before and after the boundary surface (for example, a thin film structure) of the fluid and the structure, and the flow may not be approximated accurately by the basis function of the general finite element method if the mesh mismatch occurs between the fluid and the structure. This is because the basis function of the finite element method is usually defined so that the pressure and the velocity gradient within the element become continuous functions, and the discontinuity within the element may not be represented. For example, a method has been proposed in de Hart et al., to interact a blood flow and a valva aortae (aortic valve) by applying a condition {dot over (d)}=v in which the velocities of the fluid and the structure match on the boundary of the fluid and the structure, using the finite element in which the function of pressure and the velocity gradient are elementwise discontinuous. However, this proposed method does not perform the simulation with a realistic Reynolds number. In addition, this proposed method may generate an unnatural back flow (or reflux) before and after the structure due to the effects of the constraint condition.
A technique that simulates an interaction of the fluid and the structure using the method of Lagrangian multipliers by using an incompressibility condition of the fluid in a neighborhood of the boundary of the fluid and the structure as the limiting condition has been proposed in Hisada et al., for example. However, this proposed technique neglects a viscous condition at the boundary interface due to viscosity. Hence, in the fluid structure interaction simulation method that does not depend on the mesh consistency, it may be difficult to perform an accurate simulation.
As described above, the conventional fluid structure interaction simulation method using the IT type ALE method depends on the mesh consistency between the fluid and the structure, and the load of the computations on the computer is relatively large due to the relatively large amount of computations required for the mesh control. On the other hand, in the IC (Interface-Capturing) type conventional fluid structure interaction simulation method that does not depend on the mesh consistency, it may be difficult to perform an accurate simulation. Hence, the conventional fluid structure interaction simulation methods suffer from a problem in that it may be difficult to perform an accurate simulation without being dependent on the mesh consistency between the fluid and the structure.